# Loading required package
library(margins)

######################
# AGGREGATE ANALYSIS #
######################
############
# FIGURE 2 #
############
# Loading the data
load("analysis_2_aggregate_data.RData")

### Calculating the proportions of support for obstruction for each level of intensity and each partisan group ###
# Copartisan
# Severe
severe.cop.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Severe"),c("copartisan_proportion_support")]))
# Moderate
moderate.cop.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Moderate"),c("copartisan_proportion_support")]))
# Weak
weak.cop.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("copartisan_proportion_support")]))
# Combining into vector for plotting
copartisans.proportion <- c(weak.cop.proportion, moderate.cop.proportion, severe.cop.proportion)
copartisans.percentage <- c(weak.cop.proportion, moderate.cop.proportion, severe.cop.proportion)*100

# Independents
# Severe
severe.ind.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Severe"),c("independent_proportion_support")]))
# Moderate
moderate.ind.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Moderate"),c("independent_proportion_support")]))
# Weak
weak.ind.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("independent_proportion_support")]))
# Combining into vector for plotting
independents.proportion <- c(weak.ind.proportion, moderate.ind.proportion, severe.ind.proportion)
independents.percentage <- c(weak.ind.proportion, moderate.ind.proportion, severe.ind.proportion)*100

# Outpartisans
# Severe
severe.out.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Severe"),c("outpartisan_proportion_support")]))
# Moderate
moderate.out.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Moderate"),c("outpartisan_proportion_support")]))
# Weak
weak.out.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("outpartisan_proportion_support")]))
# Combining into vector for plotting
outpartisans.proportion <- c(weak.out.proportion, moderate.out.proportion, severe.out.proportion)
outpartisans.percentage <- c(weak.out.proportion, moderate.out.proportion, severe.out.proportion)*100

### Formal t-tests of difference between strong and weak forms of obstruction for each partisan alignment group ###
# Copartisans proportion support, difference of -2.8, p= 0.5815
round((weak.cop.proportion*100)-(severe.cop.proportion*100),1)
t.test(unlist(aggregate_data[which(aggregate_data$intensity == "Severe"),c("copartisan_proportion_support")])*100,
       unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("copartisan_proportion_support")])*100)
# Moderates proportion support, difference of 19.0, p= 0.001024
round((weak.ind.proportion*100)-(severe.ind.proportion*100),1)
t.test(unlist(aggregate_data[which(aggregate_data$intensity == "Severe"),c("independent_proportion_support")])*100,
       unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("independent_proportion_support")])*100)
# Outpartisans proportion support, difference of 30.0, p= 0.002497
round((weak.out.proportion*100)-(severe.out.proportion*100),1)
t.test(unlist(aggregate_data[which(aggregate_data$intensity == "Severe"),c("outpartisan_proportion_support")]*100),
       unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("outpartisan_proportion_support")])*100)

group <- c(0.5,1.55,2.55)
par(mar=c(4.5,5,1,1), mgp=c(2,.25,0))
plot(group,independents.percentage,xlim=c(0,3),ylim=c(0,100), tck=0, cex.main=1, main="", cex.lab=1.2, cex.axis=1, 
     pch=24, cex=1.5, ylab="Average Percentage of Respondents Who Support Obstruction", xaxt="n", yaxt="n", xlab="Intensity of Obstruction", type="b", lwd=2.75)
axis(2,at=seq(0,100,by=10), cex.axis=1, 
     labels=c("0","10","20","30","40","50","60","70","80","90","100"), tck=0.01, las=1)
axis(1, at=c(0.5,1.5,2.5),cex.axis=0.9,labels=c("Weak","Moderate","Severe"))
lines(x=c(0.5,1.5,2.5), y=copartisans.percentage, pch=16, type="b", lty=2, lwd=2.75, cex=1.5)
lines(x=c(0.5,1.5,2.5), y=outpartisans.percentage, pch=15, col="grey60", type="b", lty=3, lwd=2.75, cex=1.5)
text(x=2, y=80, "Outpartisans (30.0, p<0.01)", cex=1)
text(x=2, y=16, "Copartisans (-2.8, p<0.59)", cex=1)
text(x=1.5, y=56, "Independents (19.0, p<0.01)", cex=1)

######################
# IN-TEXT DISCUSSION #
######################
### Pooling moderate and severe forms of obstruction (footnote 32) ###
# Copartisans
severe.moderate.cop.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity %in% c("Severe","Moderate")),c("copartisan_proportion_support")]))
weak.cop.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("copartisan_proportion_support")]))

# Independents
severe.moderate.ind.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity %in% c("Severe","Moderate")),c("independent_proportion_support")]))
weak.ind.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("independent_proportion_support")]))

# Outpartisans
severe.moderate.out.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity %in% c("Severe","Moderate")),c("outpartisan_proportion_support")]))
weak.out.proportion <- mean(unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("outpartisan_proportion_support")]))

# Formal t-tests of difference between strong and weak forms of obstruction for each partisan alignment group
# Copartisans proportion support, difference of -2.3, p= 0.6312
round((weak.cop.proportion*100)-(severe.moderate.cop.proportion*100),1)
t.test(unlist(aggregate_data[which(aggregate_data$intensity %in% c("Severe","Moderate")),c("copartisan_proportion_support")])*100,
       unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("copartisan_proportion_support")])*100)
# Moderates proportion support, difference of 15.0, p= 0.0002067
round((weak.ind.proportion*100)-(severe.moderate.ind.proportion*100),1)
t.test(unlist(aggregate_data[which(aggregate_data$intensity %in% c("Severe","Moderate")),c("independent_proportion_support")])*100,
       unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("independent_proportion_support")])*100)
# Outpartisans proportion support, difference of 23.0, p= 4.013e-06
round((weak.out.proportion*100)-(severe.moderate.out.proportion*100),1)
t.test(unlist(aggregate_data[which(aggregate_data$intensity %in% c("Severe","Moderate")),c("outpartisan_proportion_support")]*100),
       unlist(aggregate_data[which(aggregate_data$intensity == "Weak"),c("outpartisan_proportion_support")])*100)

### Average support for obstruction by partisan group (footnote 33) ###
# Copartisans
round(mean(unlist(aggregate_data[,c("copartisan_proportion_support")])),3)

# Independents
round(mean(unlist(aggregate_data[,c("independent_proportion_support")])),3)

# Outpartisans
round(mean(unlist(aggregate_data[,c("outpartisan_proportion_support")])),3)

#############################
# INDIVIDUAL-LEVEL ANALYSIS #
#############################
rm(list=ls())

# Loading the dataset of the HH and MC polls
load("analysis_2_individual_surveys.RData")

# Loading the dataset of all the individual surveys from iPoll
# You can generate this dataset with the parse_ipoll_data.R file
# To generate the dataset, you will need access to the Roper iPoll database
# See the instructions in the read_me_jlc_obstruction.txt file or the parse_ipoll_data.R file for more information
load("ipoll_individual_surveys.RData")

# Merging the two datasets together to create the full dataset for the analysis
omnibus_dataset <- rbind(individual_dataset, ipoll_individual_dataset)
rm(individual_dataset, ipoll_individual_dataset)

###########
# TABLE 2 #
###########
### Running the regression ###
# Model 3 of Table D.3 (used to generate predicted probabilities): gender and race covariates (available for all surveys), interaction, president fixed effects
omnibus.regression.3 <- lm(support_delay ~ copartisan_of_president*intensity + outpartisan_of_president*intensity + white + male + president, data=omnibus_dataset)
summary(omnibus.regression.3)

### Generating predicted probabilities from model 3 for the table ###
# Independents, copartisans, outpartisans for each of weak/moderate/severe intensities
# Holding other covariates constant
predict.dataset <- data.frame(copartisan_of_president=c(0,1,0,0,1,0,0,1,0), 
                              outpartisan_of_president=c(0,0,1,0,0,1,0,0,1),
                              intensity=c("weak","weak","weak","moderate","moderate","moderate","severe","severe","severe"),
                              white=c(1,1,1,1,1,1,1,1,1),
                              male=c(1,1,1,1,1,1,1,1,1),
                              president=c("bush","bush","bush","bush","bush","bush","bush","bush","bush"))
predictions  <- predict(omnibus.regression.3, predict.dataset, type=c("response"), se.fit=T)

# Outpartisan predictions for weak/moderate/severe
outpartisans <- c(predictions$fit[3],predictions$fit[6],predictions$fit[9])
round(outpartisans,3)
# Weak/severe difference
round(predictions$fit[9]-predictions$fit[3],3)

# Copartisan predictions for weak/moderate/severe
copartisans <- c(predictions$fit[2],predictions$fit[5],predictions$fit[8])
round(copartisans,3)
# Weak/severe difference
round(predictions$fit[8]-predictions$fit[2],3)

# Independent predictions for weak/moderate/severe
independents <- c(predictions$fit[1],predictions$fit[4],predictions$fit[7])
round(independents,3)
# Weak/severe difference
round(predictions$fit[7]-predictions$fit[1],3)

### Marginal effect of intensity (with weak as baseline) for each of independent, outpartisan, and copartisan groups ###
margins.omnibus.fes <- margins(omnibus.regression.3, variables='intensity', at=list(copartisan_of_president=c(0,1), outpartisan_of_president=c(0,1)))
# Severe/weak difference, for independents
summary(margins.omnibus.fes)[5,]
# Severe/weak difference, for outpartisans
summary(margins.omnibus.fes)[6,]
# Severe/weak difference, for copartisans
summary(margins.omnibus.fes)[7,]
